clc;clear;

%% 参数设置
% 固定参数（论文中的设定）
b = 0.7;         % E_j / V0
r = 0.8;         % R_j / a
c = 0.1;         % a^2 * C_j / L_j
n = 0.15;        % η
A = 2;         % 压电信号振幅

% 耦合强度
lambdaL = 10;    % 电感耦合
lambdaR = 10;    % 电阻耦合
lambdaC = 0.2;   % 电容耦合

% 外部信号频率
% 低频区间：0.002–0.07 Hz
% 中频区间：0.07–0.15 Hz
% 高频区间：0.15–0.2 Hz
% f_list = [0.002, 0.07, 0.15, 0.2];  % 单位 Hz
f_list = [0.07, 0.15];   % 只测试2个频率
% 混沌电流强度
% M_values = 0:0.05:1;
M_values = [0, 0.3, 0.6]; % 只测试3个混沌强度

% Lorenz 系统参数
sigma = 10; rho = 28; beta = 8/3;

% 时间设置
dt = 0.005;                 % 仿真步长
% transient_T = 5500;         % 瞬态时间（论文中设定为5500时间单位）
% record_T = 1000;            % 数据采集时间
transient_T = 500;   % 瞬态减少
record_T = 300;      % 采集时间也缩短
fs = 1/dt;                  % 采样频率

% 检查频率合法性：确保 f 不超过奈奎斯特频率 fs/2
if any(f_list > fs/2)
    error('存在频率超过奈奎斯特频率 fs/2 = %.2f Hz', fs/2);
end

steps_transient = round(transient_T/dt);
steps_record = round(record_T/dt);

%% 统一初始化初始状态
initial_state = [0.2; 0.1; 0.2; 0.1; 0];
init_inductive = initial_state;      % 电感耦合状态向量：[x1; y1; x2; y2; z]
init_resistive = initial_state(1:4);    % 电阻耦合状态向量：[x1; y1; x2; y2]
init_capacitive = initial_state(1:4);    % 电容耦合状态向量：[x1; y1; x2; y2]
init_lorenz = [20; 0; 0];                            % Lorenz 系统初始状态 [u; v; w]

%% 预分配存储 Q 值的矩阵
numF = length(f_list);
numM = length(M_values);
Q_inductive = zeros(numF, numM);
Q_resistive = zeros(numF, numM);
Q_capacitive = zeros(numF, numM);

%% 进度显示
total_iters = numF * numM;
progress = waitbar(0, 'Processing Fig4 simulation...');

%% 主循环：对每个外部信号频率 f 和混沌电流强度 M 进行仿真
for fi = 1:numF
    f = f_list(fi);
    for mi = 1:numM
        M = M_values(mi);
        
        % 复制初始状态（确保每次仿真条件一致）
        state_ind = init_inductive;
        state_res = init_resistive;
        state_cap = init_capacitive;
        lorenz_state = init_lorenz;
        t = 0;
        
        %% 瞬态仿真
        for k = 1:steps_transient
            % 更新 Lorenz 系统
            lorenz_state = rk4_lorenz(lorenz_state, dt, sigma, rho, beta);
            u_val = lorenz_state(1);
            % 外部声信号
            u_PC = A * cos(2*pi*f*t);
            % 更新各耦合模型
            state_ind = rk4_inductive(state_ind, dt, b, r, c, n, lambdaL, M*u_val, u_PC);
            state_res = rk4_resistive(state_res, dt, b, r, c, n, lambdaR, M*u_val, u_PC);
            state_cap = rk4_capacitive(state_cap, dt, b, r, c, n, lambdaC, M*u_val, u_PC);
            t = t + dt;
        end
        
        %% 数据采集阶段
        % 预分配存储 采集稳态信号
        x_series_ind = zeros(steps_record, 1);
        x_series_res = zeros(steps_record, 1);
        x_series_cap = zeros(steps_record, 1);
        
        for k = 1:steps_record
            % 更新 Lorenz 系统
            lorenz_state = rk4_lorenz(lorenz_state, dt, sigma, rho, beta);
            u_val = lorenz_state(1);
            u_PC = A * cos(2*pi*f*t);
            % 更新各耦合模型
            state_ind = rk4_inductive(state_ind, dt, b, r, c, n, lambdaL, M*u_val, u_PC);
            state_res = rk4_resistive(state_res, dt, b, r, c, n, lambdaR, M*u_val, u_PC);
            state_cap = rk4_capacitive(state_cap, dt, b, r, c, n, lambdaC, M*u_val, u_PC);
            % 记录中央神经元膜电位
            x_series_ind(k) = state_ind(1);
            x_series_res(k) = state_res(1);
            x_series_cap(k) = state_cap(1);
            t = t + dt;
        end
        
        %% 计算傅里叶系数 Q（调用 fourier.m）
        Q_ind = fourier(x_series_ind, fs, f);
        Q_res = fourier(x_series_res, fs, f);
        Q_cap = fourier(x_series_cap, fs, f);
        
        Q_inductive(fi, mi) = Q_ind;
        Q_resistive(fi, mi) = Q_res;
        Q_capacitive(fi, mi) = Q_cap;
        
        % 更新进度条
        curr_iter = (fi-1)*numM + mi;
        waitbar(curr_iter/total_iters, progress);
    end
end
close(progress);

%% 绘图
figure('Position', [100, 100, 1200, 400]);
colors = lines(numF);

% 电感耦合子图
subplot(1,3,1); hold on; grid on;
for fi = 1:numF
    plot(M_values, Q_inductive(fi,:), 'Color', colors(fi,:), 'LineWidth', 1.5, ...
        'DisplayName', sprintf('f=%.3f', f_list(fi)));
end
title('(a) 电感耦合 (\lambda_L = 10)');
xlabel('混沌电流强度 M'); ylabel('傅里叶系数 Q');
legend('Location', 'northeastoutside');

% 电阻耦合子图
subplot(1,3,2); hold on; grid on;
for fi = 1:numF
    plot(M_values, Q_resistive(fi,:), 'Color', colors(fi,:), 'LineWidth', 1.5, ...
        'DisplayName', sprintf('f=%.3f', f_list(fi)));
end
title('(b) 电阻耦合 (\lambda_R = 10)');
xlabel('混沌电流强度 M'); ylabel('傅里叶系数 Q');
legend('Location', 'northeastoutside');

% 电容耦合子图
subplot(1,3,3); hold on; grid on;
for fi = 1:numF
    plot(M_values, Q_capacitive(fi,:), 'Color', colors(fi,:), 'LineWidth', 1.5, ...
        'DisplayName', sprintf('f=%.3f', f_list(fi)));
end
title('(c) 电容耦合 (\lambda_C = 0.2)');
xlabel('混沌电流强度 M'); ylabel('傅里叶系数 Q');
legend('Location', 'northeastoutside');